METHOD AND APPARATUS FOR MODELING THREE-DIMENSIONAL 
ELECTROMAGNETIC SCATTERING FROM ARBITRARILY SHAPED THREE- 
DIMENSIONAL OBJECTS 
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FIELD OF THE INVENTION 
The present invention relates generally to computer aided design and object modeling. 
More specifically, it relates to a computer program product, apparatus, and method for 
detennining equivalent surface sources (or actual surface currents in the case of a metal 
20 object) on the surface of a closed, arbitrarily shaped object (e.g., a component of a 
vehicle such as an airplane part) that result from the imposition of an arbitrary time- 
harmonic incident field (e.g., a radar wave). 
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BACKGROUND 

Software utilizing boundary integral equation techniques to solve electromagnetic 
scattering problems (e.g., using the method of moments) has existed at least since the 
1960s. Such techniques have the potential for making very accurate radar cross-section 
5 predictions. However, practical application of the techniques to realistic problems is 
severely limited by poor scalability of computational resource requirements (i.e., CPU 
time and memory) relative to the problem size and complexity. With the relatively recent 
introduction of computational techniques involving high-order and "fast" mathematical 
methods (such as the fast multipole method (FMM)), the situation has changed 

. 5 s ?! 

!*( 10 drastically. Now, instead of the resource requirements for a problem described by N 

Tins' 

K unknowns scaling in 0(N 3 ) time, they scale in 0(NlogN) time, much like the Fast Fourier 

y, Transform (FFT) algorithm. Consequently, it has become possible to handle a much 

n wider array of problems of practical interest. 

if! 

;2 15 Although the methods used have improved drastically, fast methods must still be used in 
^ conjunction with iterative solvers. If the set of linear equations that arises from 

discretizing a boundary integral equation is ill conditioned, the number of iterations 
required to obtain an accurate solution is not well controlled, and may be very large. The 
overall efficiency of the solution process suffers accordingly. In addition, a high degree 
20 of numerical precision is generally required in the solution process, causing the need for 
rounding, and resulting in a loss of solution accuracy. 
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Unfortunately, the most widely used integral equation formulation for closed, perfect 
electrical conductors (PEC), namely, the combined field integral (CFIE) equation, is not, 
in general, well conditioned. The ill conditioning stems from the fact that the CFIE 
necessarily includes an integral operator component that is "hypersingular". This 
5 component is, in fact, not simply an integral operator at all, but rather the gradient of an 
integral operator. Whereas, ordinary integral operators tend to filter out the numerical 
"noise" associated with a discretization, differential operators tend to amplify it. 

The eigenvalue spectrum of a hypersingular integral operator is unbounded, which leads 
directly to an ill-conditioned linear system. The ill effects become evident in the form of 
an increase in the required number of iterations as one refines the discretization of a 
target. Fine discretizations are often required to spatially resolve source variations or 
geometric detail on the subwavelength scale. This difficulty with the CFIE (as well as 
with some other boundary integral formulations in computational electromagnetics) is 
known as the "low frequency problem/' 

Various methods have been employed to overcome this problem. The most commonly 
practiced method is to multiply the CFIE matrix by a "preconditioner" matrix at every 
step of the iterative solver procedure. An effective preconditioner improves the condition 
20 number of a linear system, thereby reducing the number of iterations required for 

solution. Examples of preconditioners that have been utilized are the block-diagonal 
preconditioner, the incomplete LU preconditioner, and the sparse approximate inverse 
preconditioner. Preconditioners such as these are general purpose - in other words, they 
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can be applied with varying degrees of success to any matrix. Unfortunately, however, 
none of them takes specific advantage of the analytical properties of the CFIE, so none of 
them is ideally suited for use with the CFIE. 

5 Instead of preconditioning the equation, it is also possible to "post-condition" it. This 
involves defining the unknown part of the original linear equation to be the driving term 
for some new linear system. If the combined linear system is properly conditioned, then 
the problem is effectively solved because, once the new auxiliary unknowns have been 
solved for, it is straightforward to determine the solution to the original linear system. 
Jjf 10 The moment method approaches that use "loop-star" basis functions fall into this 
S category. While some success has been achieved with this technique, it is cumbersome 
2 with limited extensibility, and has been applied only to low-order, moment method 
a: discretizations. 

Wj 1 5 Therefore, it is desirable to provide a preconditioner method that is specifically tailored 
*** for use with the CFIE. More specifically, the preconditioner should be designed to 

minimize the number of iterations required for a solution. This is especially desirable 
when the use of a fine discretization to resolve source variations or geometric detail on a 
subwavelength scale results in an ill-conditioned linear equation (the "low-frequency" 
20 problem). Most specifically, it is desirable to provide a computer program product, an 
apparatus, and a method for modeling the electromagnetic response of a closed arbitrarily 
shaped three-dimensional object to an arbitrary time-harmonic incident field by means of 
a well-conditioned boundary integral equation (BIE) by discretizing the well-conditioned 
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BIE to provide and solve a well-conditioned finite dimensional linear system to 
determine the distribution of equivalent surface sources on the object. 
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SUMMARY OF THE INVENTION 
The present invention provides an improved method for modeling electromagnetic 
scattering from an arbitrarily shaped three-dimensional object to an arbitrary time- 
harmonic incident field. It allows for the determination of the scattering with an 

h § 10 improved degree of numerical accuracy, while also providing for increased 

^ computational efficiency. 

Si In a first embodiment, the present invention provides a computer program product, used 

(H in a computer system having a processor, for modeling an electromagnetic response of an 
W1 15 arbitrarily shaped three-dimensional object to an arbitrary time-harmonic incident field 
by means of a well-conditioned boundary integral equation (BIE). The computer 
program product comprises recording medium; mans, recorded on the recording medium 
for producing a discretized representation of a well-conditioned BIE to provide a well- 
conditioned, finite-dimensional linear system whose solution approximates a distribution 
20 of equivalent surface sources on a model caused by the arbitrary time-harmonic incident 
field; and means, recorded on the recording medium, for solving the linear system to 
determine values of the equivalent surface sources; whereby solution of the well- 
conditioned finite-dimensional linear system provides improved computational efficiency 
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and solution accuracy. Additionally, this embodiment of the present invention may 
comprise means, recorded on the recording medium for modeling electromagnetic 
scattering from the three-dimensional object due to the time-harmonic incident field as 
determined from the equivalent surface sources. Preferably, the computer program 
5 product is tailored for determining the electromagnetic scattering from a time-harmonic 
incident field in the form of a radar wavefront and where the three-dimensional object 
represents a component of a vehicle. 

In another, more specific embodiment, the BIE produced by the computer program 
€l 1 0 product is in the form of: 

(E excitation term)+ y(H excitation term) = {(E integral operator term)+ y(R integral operation term)}j 


S| where, y is a complex constant and J represents the unknown surface source distribution. 

HI 1 5 In yet another, still more specific embodiment, the BIE produced by the computer 
C program product is in the form of: 

kT s (ik)nxE mc +iZT\ik)n-n mc = {i [-T*(ik)oT 7 \k) + T\ik)oT L (k)-k 2 T\ik)oT\k)]-K + (k)}j 

In a further embodiment, the present invention provides an apparatus for modeling the 

20 electromagnetic response of an arbitrarily shaped three-dimensional object to an arbitrary 

time-harmonic incident field by means of a well-conditioned boundary integral equation 

(BIE). The apparatus comprises a computer system having a processor, a memory 

connected with the processor, and an input/output device connected with the processor 

for receiving parameters representing an arbitrarily shaped three-dimensional object and 
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an arbitrary time-harmonic incident field; means, operating in the processor and the 
memory of the computer for utilizing the parameters representing the arbitrarily shaped 
three-dimensional object and the arbitrary time-harmonic incident field for producing a 
discretized representation of a well-conditioned BIE to provide a well-conditioned finite 
5 dimensional linear system whose solution approximates a distribution of equivalent 
surface sources on a model caused by the time-harmonic incident field; and means, 
operating in the processor and the memory of the computer for solving the linear system 
to determine values of the equivalent surface sources; whereby solution of the well- 
conditioned finite dimensional linear system provides improved computational efficiency 
10 and solution accuracy. Similar to the computer program product embodiments, the 
Jjf apparatus further comprises means, operating in the processor and the memory of the 
12 computer for modeling electromagnetic scattering from the three-dimensional object due 
S! to the time-harmonic incident field as determined from the equivalent surface sources. 

m Again, the apparatus is preferably tailored for determining the electromagnetic scattering 
Wl 1 5 from a time-harmonic incident field represented in the form of a radar wavefront and 
|s * where the three-dimensional object represents a component of a vehicle. 

In another, more specific embodiment, the BIE produced by the apparatus is in the form 
of: 

20 (E excitation term)+ /(H excitation term) = {(E integral operator term)+ y(H integral operation term)}j 
where, y is a complex constant and J represents the unknown surface source distribution. 


Page 9 of 55 
6/15/2001 
HRL053FINAL 


In yet another, still more specific embodiment, the BIE produced by the apparatus 
product is in the form of: 

kT s (ik)n x E mc + iZr (ik)n • H mc - jf [- T a (ik) o T T (Jt) + T p (ik) o T L (k) - k 2 T s (ik) o T s (k)]~ K + (k)}j 


The present invention further provides a method for modeling an electromagnetic 
response of an arbitrarily shaped three-dimensional object to an arbitrary time-harmonic 
incident field by means of a well-conditioned boundary integral equation (BIE). The 
method operates on a computer system having a processor, a memory connected with the 
processor, and an input/output device connected with the processor for receiving 
parameters representing an arbitrarily shaped three-dimensional object and an arbitrary 
time-harmonic incident field. The method comprises the steps of: utilizing the 
parameters representing the arbitrarily shaped three-dimensional object and the arbitrary 
time-harmonic incident field for producing a discretized representation of a well- 
conditioned BIE to provide a well-conditioned finite dimensional linear system whose 
solution approximates a distribution of equivalent surface sources on a model caused by 
the arbitrary time-harmonic incident field; and solving the linear system to determine 
values of the equivalent surface sources; whereby solution of the well-conditioned finite- 
dimensional linear system provides improved computational efficiency and solution 
accuracy. 

Additionally, the method can further comprise the step of modeling electromagnetic 
scattering from the three-dimensional object due to the time-harmonic incident field as 
determined from the equivalent surface sources. The method is preferably operated 
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where the time-harmonic incident field is a radar wavefront and wherein the three- 
dimensional object represents a component of a vehicle, in order to determine the 
scattering of radar waves from a vehicle component due to in incident radar wave. 

5 In another, more specific embodiment, the BIE produced by the method is in the form of: 

(E excitation term)+^(H excitation term) = {(E integral operator term)+ y(R integral operation term)} J 

where, y is a complex constant and J represents the unknown surface source distribution. 

% 10 In yet another, still more specific embodiment, the BIE produced by the method is in the 
jjg form of: 

j| kT s (ik)n xE ,K + iZT\ik)n • H'" 0 = jf [- T*(ik) ° T r (k) + T f (ik) ° T L (k) - k 2 T s (ik) ° T s (k)]- K + (k)}j 

j| DESCRIPTION OF THE FIGURES 

in 15 FIG. 1 is a block diagram representing a general purpose computer system used in 

Is? ?. 

! J conjunction with the present invention; 

FIG. 2(a) is an illustrative diagram showing an object having a diameter of 2.5X (where X 
represents a wavelength), and corresponding to element 300 of FIG. 3; 
FIG. 2(b) is an illustrative diagram showing the object of FIG. 2(a) modeled as a plurality 
20 of triangular patches, and corresponding to element 302 of FIG. 3; 

FIG. 2(c) is an illustrative diagram showing the object of FIG. 2(b) with an incident field, 
and corresponding to element 304 of FIG. 3; 
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FIG. 2(d) is an illustrative diagram showing the object of FIG, 2(c) with an incident field 
and with electromagnetic fields resulting from equivalent surface sources due to the 
incident field, and corresponding to elements 306, 308, and 310 of FIG. 3; 
FIG. 3 is a flow diagram demonstrating the method of the present invention in 
5 conjunction with the steps of inputting object and incident field parameters, dividing the 
model of the object into a plurality of patches, and the optional step of determining 
electromagnetic fields; and 

FIG. 4 is an illustrative diagram of a computer program product according to the present 
invention, stored on CD or DVD. 

i 10 

| DETAILED DESCRIPTION 

y* The present invention provides a computer program product, a method, and an apparatus 

"a. S 

s . m for generating and solving well-conditioned boundary integral equations for modeling 
if! three dimensional electromagnetic scattering of an incident wavefront from an object. It 
W 1 5 should be understood that this description is not intended to limit the invention. On the 
contrary, the invention is intended to cover alternatives, modifications and equivalents, 
which are included in the spirit and scope of the invention as defined by the appended 
claims. Furthermore, in the following detailed description of the present invention 
numerous specific details are set forth in order to provide a thorough understanding of the 
20 present invention. However, it will be obvious to one of ordinary skill in the art that the 
present invention may be practiced without the specific details. 


Introduction 


Page 12 of 55 

6/15/2001 
HRL053FINAL 


In the presence of an incident wavefront, the surface properties of an object affect 
scattering of the wave from the object's surface. By properly modeling the object and the 
incoming wavefront, it is possible to determine how waves will be scattered from the 
object's surface. Determination of this scattering is of particular importance with regard 
5 to designing objects having reduced radar cross-sections as well as in designing antennae 
or other objects that are configured to maximize electromagnetic interception. 

The present invention, in a general sense, provides a computer program product, a 
method, and an apparatus for receiving a computer model of an object comprised of a 
*Q 1 0 mesh of a plurality of patches; for specifying an incident wavefront having an excitation 
W wavelength, wavefront geometry (e.g., plane wave, etc.), and incident angle; and, based 
on the properties of the object and the wavefront, determining the equivalent surface 
sources resulting from the interaction of the wavefront with the object, and from the 
S equivalent surface sources, modeling the scattering of energy from the surface of the 

m 15 object 
sub- 
invention Overview 

The following overview is presented with reference to the figures in order to provide a 
high-level conceptual understanding of the present invention prior to discussing the 
20 mathematical details and specific advantages inherent therein. First, an overview of the 
data processing system of the present invention is provided, followed by a brief 
illustration and flow diagram explaining the operation of the present invention. 
Subsequent to the overview, further details of the present invention will be provided, 
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followed by a set of numerical results demonstrating the advantages of the modified 
CFIE of the present invention over the conventional CFIE. Finally, an appendix and 
glossary are provided to give the reader a convenient source for understanding some of 
the terms and conventions utilized herein. 

5 

Data Processing System 
The present invention is designed to be operated on a data processing system, which may 
comprise a general or special purpose computer. Having reference now to the drawings, 
in FIG. 1, there is shown a computer or data processing system generally designated by 
m \0 reference character 100. The computer 100 generally includes a central processing unit 
!ff (CPU) 102 which is coupled with a memory 104. The memory 104 generally comprises 
}*? both random access memory (RAM) and read-only memory (ROM) portions. The 

computer 100 further comprises at least one input/output device 106 connected with the 
|p CPU 102, such as a keyboard, a pointing device, a scanning device, or other available and 
III 1 5 useful devices. A display 108 is also connected with the CPU 102 to provide output to a 
\*k user Further a storage mechanism 110 is provided to enable a user to store and retrieve 
data. Non-limiting examples of various storage mechanisms 110 that may be of use with 
the present invention include hard drives, optical drives (such as a compact disc or digital 
versatile disc), and magnetic drives (such as floppy disks, magnetic tapes, and Zip™ 
20 disks). Finally, a communications adapter 112 (such as a modem or network interface 
card) may also coupled with the computer 100 in order to allow communication with 
other computers 100. Note further that the data processing system may also comprise a 
network of computers either local to or remote from each other. In addition, the various 
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components of the data processing system may all be local or may be distributed 
remotely from each other. 


Although described at a major-component level, it is assumed that the computer 100 
5 includes all of the necessary adapters to drive all of the peripherals necessary for use with 
the present invention. Additionally, various commercially available processors could be 
used in conjunction with the present invention, and could be programmed to operate the 
logic of the preferred embodiments. 

\Q 1 0 Illustration and Flow Diagram of the Operation of the Present Invention 

ffl Referring to FIG. 2, and FIG. 3, the operation of the present invention is described by 
means of an illustrative description and a flowchart. After starting operation of the 
present invention, represented by the begin 300 element in FIG, 3, an object 200 is 
5 divided into a plurality of patches, represented by element 202 in FIG. 2 and element 302 
U115 in FIG. 3. Note that methods for dividing the object into a plurality of patches are well- 
H ; known in the literature, and that the patches may be of an arbitrary geometric shape (e.g., 
triangles, squares, etc.). Next, parameters are supplied for the incident field and the 
properties of the object as depicted by element 304 of FIG. 3. The incident field is 
depicted as element 204 in FIG. 2. The dashed lines surrounding elements 302 and 304 
20 of FIG. 3 are intended to indicate that they are environmental to the present invention, as 
the patches 202 representing the object 200 and the parameters for the incident field 204 
and the object 200 may simply be passed as inputs to a system incorporating the present 
invention. Optionally, however, the present invention may include a system for dividing 

Page 15 of 55 

6/15/2001 
HRL053FINAL 


the object 200 into a plurality of patches 202 and for more directly determining the 
associated parameters. In addition, although the object 200 is represented as a sphere and 
the incident field 204 is represented as a series of planar waves, these specifics are 
provided for illustrative convenience only. The object 200 may be any closed arbitrary 
5 shape and the incident field 204 may be any arbitrary time-harmonic field. Next, as 
shown by element 306 in FIG. 3, the system develops a discretized representation of a 
well-conditioned boundary integral equation (which takes the form of a well-conditioned 
matrix equation), modeling the interaction between the incident field 204 and the object 
200, The well-conditioned linear system is then solved numerically to determine 

<w 

iiQ 10 equivalent surface sources on the object, as depicted by element 308 in FIG. 3. After 

13 solving for the equivalent sources, the present invention may utilize the equivalent 

W sources in order to determine electromagnetic fields resulting therefrom, as depicted by 

^ element 310 of FIG. 3. The dash-dot-dot line in element 310 is intended to indicate that 

m it is optional. The electromagnetic scattering from the object is represented by element 

m 15 206 in FIG. 2. 

U 

By providing a well-conditioned linear system in the determination of the equivalent 
sources, the present invention helps to ensure that the solution to the linear system can be 
performed with fewer iterations and with increased numerical precision relative to 
20 conventional methods. In fact, without the use of the well-conditioned linear system, it is 
sometimes impractical to determine the equivalent surface sources accurately due to the 
limited numerical precision of one's computer. 
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Recording Medium 

In one embodiment, the present invention takes the form of a computer program product 
recorded on a recording medium. A disc, shown in FIG. 4 represents a recording 
medium 400 such as a compact disc (CD) or a digital versatile disc (DVD). Although 
5 depicted as a CD or DVD, any recording medium could be utilized, examples of which 
include magnetic discs and tapes, and other optical storage medium. 

Further Details of the Present Invention 

Well-Conditioned CFIE Formulation 
*f| 10 In accordance with features of the present invention, to provide improved computational 

jfTT; 

m efficiency and accuracy, modifications are made to the combined field integral equation 
:^ (CFIE) for the solution of electromagnetic scattering problems in three dimensions. 

" Unlike the conventional CFIE, the modified CFIE of the present invention is well- 
m conditioned. It combines, a standard magnetic field integral operator with a pre- 
fH 15 conditioned electric field integral operator, with the resulting modified CFIE being an 
M integral equation of the second-kind that has no spurious resonances. As a result of the 
pre-conditioning, the number of iterations required to solve the CFIE stabilizes even as 
discretization of the surface is refined, whereas refinement of the discretization in 
conjunction with the conventional CFIE results in an ill-conditioned linear system and a 
20 growing number of iterations. The conventional CFIE, not being well-conditioned, 
suffers from two major drawbacks. 
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As mentioned previously, there are three principal reasons why well-conditioned 
boundary integral formulations are beneficial The first is the fact that sparse 
decompositions of the integral operators of scattering theory, but not their inverses, exist, 
so iterative solvers must be employed. Well-conditioned systems of equations can be 
5 solved with few iterations. Second, using a fine discretization to resolve source 
variations or geometric detail on a sub-wavelength scale typically results in an ill- 
conditioned linear equation as well (also known as the "low-frequency" problem in 
computational electromagnetics). Third, because the process of generating a solution to 
the set of discretized equations can require a very wide range of values (in terms of 

IKS? 

?g 10 magnitude), the resulting solutions may be quite inaccurate because a high degree of 
W numerical precision may be required, which may exceed the accuracy permitted by the 

fil. 

data processing system. 

|jt Only second-kind integral equations, or objects with similar spectral behavior (such as 
lit 1 5 appropriately preconditioned differential equations) can be solved with fully controlled 
M approximation error. The correct operators are the sum of a constant (or at least well- 

conditioned and easily invertible) operator and a compact operator. 

Boundary integral operators of scattering typically violate this requirement in one of three 
20 ways: the spectrum may accumulate at zero; the operator may have an unbounded 

spectrum; or the operator may have small eigenvalues associated with resonances, often 
unphysical (e.g., "spurious resonances"). 
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For electromagnetic scattering from perfectly conducting (PEC) surfaces, the standard 
boundary integral equations are the electric field integral equation (EFIE) 


■ n x E mc = Ti (where inc stands for incident), (1) 


and the magnetic field integral equation (MFIE) 


ZnxET = 


1 + * 

2 J 


J, (2) 


IS 10 where the integral operators T and K are defined by 
m TJ = T(k)J 

« m 3 ikn(x) x jds'^G(k; x, x') J(x) + ~ V[VG(t, x, x') • J(x')]j and (3) 

KJ = K(k)J = - n(x) x jds'VG(k; x, x') x J(x') , (4) 


where V denotes differentiation with respect to x, and n(x) is the unit normal to the 
15 surface at x. Note that the other terms follow the usual conventions: JsZnxH is the 
unknown surface current, E'" c and YL mc are the incident electric and magnetic fields, 

respectively, Z = ^jju/s is the wave impedance, and G(k; x, x') = exp(ikr) 1 4xr is the 3d 
Helmholtz kernel with r = |x - xj being the distance separating field and source points. 
Harmonic time dependence e is assumed. 

20 


Page 19 of 55 

6/15/2001 
HRL053FINAL 


The MFIE is a second-kind integral equation. Unfortunately, this equation is suitable for 
an unacceptably small class of electromagnetic problems. It is inapplicable to open 
surfaces, becomes ill-conditioned in the presence of geometric singularities, and suffers 
from spurious resonances. The EFIE has both a compact piece and a hypersingular piece 
5 (coming from the double gradient term). The spurious resonances of the MFIE may be 
eliminated by adding the EFIE to form a combined field integral equation (CFIE). The 
cost of doing so is the introduction of EFIE 5 s hypersingular piece, which spoils the 
conditioning for fine discretizations (or low frequencies). 

10 It is known that a hypersingular integral operator and a first-kind integral operator are 

m ideal preconditioned for each other, in the sense that the composition of the two has the 

\ y 

\*j spectral characteristics of a second-kind integral operator. A similar approach may be 

" H employed to analytically precondition the EIFE. In fact, as disclosed herein, the electric 

j5 field integral operator T preconditions itself 

in 15 

|U Two issues are important for the successful application of this idea to closed bodies. 

First, only the local (or short distance) behavior of the preconditioner is important for 
asymptotic conditioning. Thus, the EFIE may be preconditioned by multiplying it by an 
electric field integral operator corresponding to an arbitrary wavenumber, real or 
20 complex; if the wavenumber has a positive imaginary part, the introduction of any 

additional resonances may be avoided. If the EFIE preconditioner reproduced the MFIE 
resonances, then the CFIE would also have them. Second, one must take care that the 
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discretization of the product of the preconditioner and preconditioned operators preserves 
the correct spectral properties. 


The present invention provides well-conditioned formulations for closed surfaces. 
5 Example numerical results are provided further below in order to demonstrate the 
advantages of the modified CFIE formulation over the conventional CFIK 


In the literature, integral operators have been constructed from the kernel for the Laplace 
and Helmholtz equations in 2 d. Through the construction of these operators, it has been 
i Jj 10 observed that the product of a first-kind operator, constructed from an undifferentiated 

ro kernel, and a hypersingular operator, constructed from a twice differentiated kernel, has 

ftf 

\& the desirable spectral characteristics of a second-kind operator. Since the EFIE integral 

^ operator T has both of these, it may be seen that the composition of two such operators 

|S t 2 = T o T would include a constant operator and a compact operator. Also of concern is 

IH 1 5 that the product of hypersingular components would produce another hypersingular 
H operator. However, it is easy to see that the rotation operation n x in the definition of T 9 

as shown in equation (3) above, annihilates the component of the surface vector field 
normal to the surface and ensures that the product of the two hypersingular operators is 
identically equal to zero. Indeed, applying the hypersingular component of the second T 
20 operator to an arbitrary tangential surface vector function f(x 5 ) produces a surface 
gradient function, 


nxVj^(x) - -(nx v) Jds' VG(k; x, x') -f(x') , (5) 

* s 
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(note the dot-product operator, •) which the hypersingular component of the first T 
operator, in turn, annihilates (for closed surfaces) by virtue of the identity 


V s -[nxV^(x)] = 0, (6) 
5 with V s denoting the surface gradient operator on S; the identity shown in equation (6), 
the surface analog of the 3d identity V • [V x (j>{x)\ = 0 , can be found and is valid for any 
sufficiently smooth function ^on S. It follows immediately from equations (3), (5), and 
(6) that T 2 behaves as a second-kind integral operator. 

IKS? 

§ 1 0 In this specification, the spectral properties of the EFIE and MFIE integral operators and 

IB 

W combinations thereof for the PEC sphere, a simple 3d target for which the spectral 
properties of these operators are known analytically, are discussed. However, this 
m discussion may be expanded to incorporate other object shapes. The use of the PEC 
W sphere herein is intended to provide a non-limiting example of an object to which the 
O 1 5 present invention may be applied. A complete set of basis functions on the surface of a 
sphere of radius a is given by the vector spherical harmonics, 

U /m (^) = nxX /m (0,p), (8) 
defined here in terms of the scalar spherical harmonics Y lm (9,<p) . 

20 The result of applying T and K + = (K + to each basis function is: 


Page 22 of 55 
6/15/2001 
HRL053FINAL 


HJ 1j; MET, (fc*)X fa J 


and 


fx, j r a;(to)H ; (^)x /ra 


,m =i v k (10) 

where J { and H, are Riccati-Bessel and (first-kind) Riccati-Hankel functions of order /, 

and k is the wavenumber associated with the kernel of each integral operator. The 
Riccati-Bessel and Riccati-Hankel functions are defined in terms of spherical Bessel and 
Hankel functions j t {x) and hf\x) by 

= x/^^and (11) 

M t (x) = xh®(x) (12) 

Although chosen basis functions Xi m and U/ OT are not eigenfunctions of the operator T{k), 
they are eigenfunctions of T 2 (k) = T(k) o T(k) : 

Tl{k ^} = " J ' (fc,) H < ( ^ } J < H ' {ka) {u*} • (13) 

The operator T 2 ^) has a bounded spectrum, since, in the limit of large /, its eigenvalues 
accumulate at -1/4 (a result which follows from the asymptotic properties of ji and h\ x) ). 

However, as is evident, from equations (10) and (13), the operator T 2 (k) also shares 
resonances (at the zeros of J\ (ka) for the Xi m modes, and at the zeros of J, (ka) and for 

the U/ m modes) with the MFIE operator K + (k) . This is also evident from the identity 
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5 


T 2 (k) = i - K 2 (k) = -K-(k)oK + (k), (14) 

where K' =--K. Therefore, although T\k) is a second-kind integral operator, it is 
2 

not a suitable component of a resonance-free combined field integral equation for closed 
bodies. 


As previously mentioned, boundedness of the spectrum of the product of two EFIE 
operators (of the form shown by equation (3)) is assured if they have the same short- 
distance behavior, a condition that does not require the two operators to share the same 

o 

5 wavenumber (propagation constant). If EFIE operators are chosen with different 
i 10 wavenumbers, T{k x ) and T(k 2 ) , it is possible to simultaneously obtain a bounded 
product and avoid MFIE resonances. 


The following analysis indicates that ik is a particularly good choice for the wavenumber 
in the preconditioning operator (assuming that the wavenumber k is real). The 
1 5 eigensystem for T (ik) ° T(k) on a sphere is 

r ^ frJM \ f > {ika) H '< {ika) J ' {ka) ^ M (15) 
T {ik) - VWp^ | = -| J/ (ika) H, (ika) 3] {ka) H; (to) U /m J ' 

It is straightforward to show (given the properties of/7 and h\ x) ) that the eigenvalues of 
T (ik) o T(k) accumulate at ^ and - l - for the X/ m and U /w eigenmodes, respectively, 
and that T (ik) o J(Jt) does not share any resonances with the MFIE operator K + (k) . 

20 
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Since T (ik) ° T{k) is a second-kind integral operator (in the sense described below) and 
does not share any resonances with K + (k) , it is possible now to develop the well- 
conditioned CFIE operator. The simplest form of such an operator is 


5 T(ik)oT(k) + aK + (k), (16) 

where a is a constant to be chosen. In developing this modified CFIE operator, the EFIE 
portion has been modified before being added to the MFIE portion (which is already a 
second-kind integral operator). The same applies to the excitation side of the equation. 

Jf 10 The resulting CFIE may be written as 

y 

y -T(ik)(nxE mc )+aZnx11'" c = [T(ik)oT(k) + aK + (k)]f . (17) 

» s The eigensystem for the CFIE operator in (1 6) is 

J l u /J I Ui ( ika ) H / i ika ) J / ( ka ) + ia J / ( ka )] H ^ ( ka ) 11 1> 


15 

If one chooses a = ±1 then, as a function of the argument ka, these eigenvalues have no 
zeros. For a = +1 , they circle the origin of the complex plane. 

Other well-conditioned (modified) CFIE operators can be devised, for example, by 
20 preconditioning the MFIE part before combining it with the preconditioned EFIE part. 
Two example forms of alternative well-conditioned (modified) CFIE operators are as 
follows: 
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T(ik)oT(k) + aK\ik)oK + (k) 
and 

T(ik) o T(k) + ccnxK + (ik) °nxK + (k). 


(19) 


(20) 


It has been observed that the numerical behavior of the three modified CFIE operators (in 
(16), (19), and (20)) is similar. 

By virtue of the analysis on spheres, the modified CFIE operators in (16), (19), and (20) 
have been shown to be second-kind and resonance-free. However, given that the 
asymptotic behavior of the eigenvalues on a smooth surface stems from the short distance 
behavior of the kernel, it may be seen that the asymptotic behavior of the various 
operators on spheres should also obtain for any closed surface that can be obtained by 
smooth deformation of a sphere. The numerical results presented further below support 
this argument. Results are also presented for a cube, which, like many targets of practical 
interest, has geometric singularities. These suggest that the modified CFIE formulations 
should be well-conditioned for a wide class of closed surfaces. 

Discretization and an Alternate Form of the Modified CFIE Formulation 

Discretization of the modified (well-conditioned) CFIE formulation utilized in the present 

invention may take place by any appropriate method. Suitable methods are common in 

the literature, non-Umiting examples of which include the Nystrom discretization method 

and the method of moments. A particularly useful reference in regard to the Nystrom 

discretization method, which may be used in conjunction with the present invention may 
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be found in reference [13], in the Background, titled "Numerical solution of the 
Helmholtz equation in 2d and 3d using a high-order Nystrom discretization" by Canino et 
al. Note also, incorporating the approach of the present invention into existing "fast" 
solvers is completely straightforward. 

5 

There are several ways to discretize the product operator T (k x )° T(k 2 ) . The simplest 
and most straightforward approach, multiplying the discretized representations of the 
individual operators, can lead to numerical difficulties. The reason is that it is relatively 
difficult to make the discretized representations of the hypersingular part of each operator 
10 sufficiently accurate (especially for high-spatial-frequency eigenmodes) to numerically 
effect the cancellation that obtains analytically. 

Effective discretizations of T (kj <=> T(k 2 ) can be obtained either by discretizing the 
product operator directly or by reformulating the product operator to eliminate the 
1 5 product of hypersingular operators. The second approach has been implemented using a 
reformulated product operator that eliminates all instances of hypersingular operators. A 
short derivation of the reformulated equation is provided below. 


The first step toward obtaining a more useful form of the product operator T (k } ) o T(k 2 ) 
20 is to separate each integral operator into its singular and hypersingular components. 
Introducing the abbreviations, 

T x =T(k 1 ) and (21) 

T 2 =T(k 2 ), (22) 
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and substituting, the following can be written: 


7; =ik{T x s +—T l 


1 » 


(23) 



iff 
2 -> 


(24) 


where, 


(25) 


7„f J = (n(x) x V) fds' VG(k m ; x, x') • J(x'). 


(26) 


5 


The product operator T x oT 2 can be expanded into four terms. Each of the two cross 
terms, T x s ° T 2 H and T X H o r/ , can be transformed (by Stokes's theorem) into a product of 
new, single-gradient integral operators on S plus a line integral around the boundary of S. 
The term formed by the product of the hypersingular integral operators, T X H ° Tf , reduces 
to a line integral. The result is further simplified by noticing that two of the three line 
integrals, when applied to J, can be combined into a single term whose argument is 
identical to the incident electric field E mc by virtue of equation (1). 

The next step is to reformulate the excitation side of the equation, taking advantage of the 
fact that the incident wave obeys Maxwell's equations. By applying Stokes's theorem, 
the term 7^[n(x 1 ) x E'" c (x')] may be rewritten as the sum of a single-gradient integral 
operator on V x E mc (x f ) and a line integral that exactly cancels the line integral 
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involving E mc on the other side of the equation. A further simplification follows from 
Faraday's Law, V x E = ico^M. . 

The final result for the analytically preconditioned EFIE with reformulated integral 
operator product is 

-/^(nxE')-Z^V(n.H') = (f ^ < + f <» " Wf • -f 2* o J?)J 

(27) 

where the various integral operators are defined by: 

T:<f> = n(x) x jds'VG(k m ; x, x'^x 1 ), (28) 

S 

T^ = n(x)x ^•n( X >V l G(i;x ) x l V(x l ) ! (29) 

T L m l = f^'VGfcx^-ftx'), (30) 

-n(x)- J^'VGfex^xftx'), (31) 

S 

Tjf = n(x) x J**G(£ m ; x, x')f (x') , and (32) 

s 

T*j = n(x) x J<a'G(*„; x, x')^(x') , and (33) 

as 

where m=\, 2. Note that T£ , T% , and T* map scalar functions to surface vector 
functions, whereas T± andrj do the reverse. The operator on the right hand side of 
equation (27) maps surface vector functions into surface vector functions. 
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Now, closed surfaces will be discussed, and it may be observed that T x o T 2 behaves like 
a second-kind integral operator. 


If S is a closed surface, the term Tf ° T 2 L J vanishes and equation (27) simplifies to 
5 -^(nxEVZ-^rCn-H') = k 2 )3, (34) 
where 

S n = S(k v k 2 ) - ° 2? + f r/ o - *^2* o T . (35) 
O Now, several features of S 12 will be discussed. 

jS 1 0 First, all of the individual operators that comprise Sn involve kernels with one or no 
!* gradients on the Helmholtz Green' s function G. All such integral operators are bounded. 

M 

M Second, the eigenvalues of the integral operator S 12 do not accumulate at the origin. This 

P may be demonstrated by examining its three components, T* o J 7 / , Tf o 7/ , and 
15 7;* o 7/ . The operator T x a ° r/ is a second-kind integral operator for the transverse 
(divergence-free) component of J, and is identically zero for the longitudinal 
(irrotational) component of J. Likewise, the operator Tf o T 2 L is a second-kind operator 
for the longitudinal component of J, and is identically zero for the transverse component 
of J. Since any surface current distribution can be decomposed into longitudinal and 

k k 

20 transverse components, the sum —T x a o T? +j~T/o T 2 L is a second-kind integral 

1 2 
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operator; subtracting k x k 2 T* ° T 2 S , a compact operator, does not change this result. As 
observed previously, resonance sharing may be avoided by setting k x = ik and £ 2 = £ . In 

this case, the eigenvalues of S 12 accumulate at two points, ± % - , rather than only at - ^ . 

Third, the spectrum ofS 12 , after discretization, is bounded and includes accumulation 
points at the expected locations. However, an accurate discretization will have zero (or 
very small) eigenvalues wherever the EFIE operator T(k 2 ) has a resonance. Thus, it has 
to be combined with an appropriate discretization of the MFIE operator, to obtain an 
effective discretization of the CFIE. 

Finally, it should be noted that equation (34) is manifestly insusceptible to the "low- 
frequency" problem that plagues the EFIE. Since the well-conditioned behavior of S 12 

comes from the composite operators ^T? °T 2 T and ^Tf <>T 2 L , both of whose 

prefactors have modulus unity (assuming |^| = \k 2 \ = k), and since the term 
k^Tf o T 2 tends to zero as k -» 0 , the full operator S n remains well conditioned in the 
limit of low frequency. 

In summary, although the operators T x ° T 2 and S I2 have identical spectral properties for 

closed bodies, it is easier to construct an accurate discretization for S }2 because it is 

20 composed of less singular integral operators. Matrix representations of S I2 have bounded 

spectra, but also suffer from spurious resonances inherited from the EFIE operator T(k 2 ) . 
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These resonances can be eliminated by combining Sn with K* (k 2 ) (or the modified 
MFIE operators in equations (19) and (20)). The result is a well-conditioned system of 
linear algebraic equations. 


5 An Alternate General Formulation of the Modified CFIE 

The boundary integral equation used for the present invention may be formulated more 
generally, as follows: 

(E excitation term)+ y(H excitation term) = {(E integral operator term)+ r(H integral operation term)}j 

O 10 

|jj where, 

m y is a complex constant, 

% J represents the unknown surface source distribution, and 

0 the (E excitation term) and (E integral operator term) are of a formulation selected from 

y 15 the group consisting of A and B, with 

;™ 

y, A representing the combination of 

(E excitation term) = T(^)(- n(x) x E'" e (x)) and 

(E integral operator term) = Tfa ) T(k 2 ) , and with 
B representing the combination of 
20 (E excitation term) = 

i^iW-nix) x E mc (x)) - (*i)(Zii(x) • H'" e (x)) and 

"i 

(E integral operator term) = 
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r^T a (t) o T T (k 2 ) + ^T\k 2 ) o T L (k 2 ) - KhT S (K) ° T s (k 2 ) ; and where 

the (H excitation term) and (H integral operation term) are of a formulation 
selected from the group consisting of A', B', and C, with 
A' representing the combination of 

(H excitation term)= Zn(x) x H'" c (x) and 

(H integral operation term) = K + (k 2 ) ; 
B' representing the combination of 

(H excitation term)= ^ + (^)(Zn(x) x H ,nc (x)) and 

(H integral operation term)= K + (*,) o K + (k 2 ) ; and 
C representing the combination of 

(H excitation term)s n x K + (k,)(Zn(x) x H'" c (x)) and 

(H integral operation term) = [nx /T (A^Jo [nx K + (k 2 )\, wherein; 
T(k) = ikn(x) x jds ' |g(*; x, x') J(x) + V[VG(£; x, x') • J(x')]J , 

.£ + (jfc) and 

K{k) = -n(x) x jds' VG(£; x, x') x J(x') , and 

s 

wherein J = Zn x H is the unknown surface current; 
£ is a complex number; 

E mc and H mc are the incident electric and magnetic fields, respectively; 

Z = «Jju/e is the wave impedance; 

Page 33 of 55 

6/15/2001 
HRL053FINAL 


G(k; x, x') = exp(ikr)/4nr is the 3d Helmholtz kernel, 

where r = |x - x'| is the distance separating field and source points with 

harmonic time dependence e' im assumed, and where 

T a {k)<t> = n(x)x jds'VG(k; x, x'^(x') , 

s 

T p (k)<f> = n(x) x J^'n(x') x VG(k; x, xW) , 
s 

T L (k)f = \dsWG(k; x,x'X(x'), 

s 

T T (k)f = n(x) • jds'VG(k; x, x') x f (x') , and 

s 

T s (k)f s n(x) x jds'G(k; x, x')f (x') , and where 

s 

V and V denote differention with respect to the field point x and source 
point x', respectively, and n(x) and n(x')are outward-pointing unit 
vectors normal to the target surface S at x and x' , respectively, and 

k 2 is the propagation constant in the medium external to S; 
k x is a complex number chosen to avoid resonances. 

This general form may be modified to provide the preferred form of the modified CFIE, 
as previously mentioned, 

/fcJ s 01)nxE'" c + /Zr a (^)n-H' KC = \i[-T\ik)oT T (k) + T\ik)oT L (k)-k 2 T s (.ik)oT s (k)]-K + (k)}j 
, where (ik) is substituted for (k) in the various parts of the equation as defined above. 

Specifically, for example, the term 
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T a {k)(f> = n(x)x fa'VGik; x, x'^(x') is simply changed to 

s 

T a (ik)<j> = n(x)x Jtfc'VGC*; x, x'V(x') . 

s 

The same type of change is applied to all of the other terms. 

Numerical Results Demonstrating the Advantages of the Modified CFIE Over the 
Conventional CFIE 

Below, the computational performance of the conventional CFIE (referred to below as 
CFIE1), 

-nx(nxE') + ZnxH' = [nx T(k) + iT (36) 
is compared with the preconditioned CFIE (referred to below as CFIE2) 
kT s {ik)n xE'+ iZT a (/Jfc)n H' -ZnxH' 

= {i [-T a (ik)oT T (k) + T^(ik)oT L (k)-k 2 T s (ik)oT s (k)\-K + (k)}j (37) 
produced by combining equation (17) (with a = -1) and equation (34) (with k x = ik and 
k 2 = k ). Individual operators in these equations were discretized using a high-order 
Nystrom scheme. In all cases, the wave impedance Z was set to unity. 

Three examples are provided below to demonstrate the operation of the present invention. 
The first example shows how the condition number of each operator, defined as the ratio 
of the largest to smallest singular values, depends on the fineness of discretization. Table 
1, below, lists the condition number (CN) of the matrix representing each CFIE operator 
as the size of the sphere decreases. In all cases, the same discretization was used, created 
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by placing a 6-point quadrature rule on each of the 80 nearly identical patches that cover 
the sphere, for a total of 960 unknowns. As the sphere radius decreases, the condition 
number for the CFIE2 integral operator stabilizes at about 2, whereas the condition 
number for the CFIE1 continues to grow in inverse proportion to the radius. 


radius(X) 

CFIE1 

CFIE2 

1 

4.2 

3.04 

1/4 

15 

2.68 

1/16 

59 

2.04 

1/64 

230 

1.99 

1/256 

940 

1.97 

1/1024 

3800 

1.97 

1/4096 

15000 

1.97 


Table 1: Condition number of CFIE matrices for shrinking PEC spheres 

The second test compares iterative solver performance for the modified CFIE with that 
for the conventional CFIE. The target geometry consists of two PEC spheres, one with a 
radius of 2 / 2 , the other set at a resonant radius, namely, the first zero of J[ (2tct I X) or 
r » 0.436674572 . The spheres are separated by a 2/1 00 gap. The patches near the gap 
were subdivided by a factor of 10 to adequately resolve the currents, which vary rapidly 
there. Table 2 compares the iteration counts and radar cross section (RCS) errors for 
several discretizations. The iterations columns list the maximum number of iterations a 
conjugate gradient squared (CGS) route required to reach a residual error of 10" 3 . A 
solution computed from a substantially more refined discretization provided an accuracy 
reference. The stated error is the root mean squared (RMS) value of the difference 
between the monostatic # RCS of the comparison solution and the reference solution at 

1 8 1 angles. For identical discretizations, the two methods had about the same error. The 
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data show a dramatic difference, however, in the iteration count behavior of the two 
methods in response to discretization refinements. 


unknowns 

patches 

CPIE1 

CPIE2 



iterations 

error 

iterations 

error 

1496 

748 

60 

0.46 

9 

0.40 

4488 

748 

126 

0.18 

11 

0.18 

996 

498 

44 

0.61 

9 

0.48 

2988 

498 

103 

0.23 

12 

0.16 

5976 

498 

163 

0.016 

11 

0.029 


Table 2: Iteration count and solution error vs, discretization for two PEC spheres. 

The third test also compared iterative solver performance for the two CFIE formulations. 
In this case, the target was a cube of size \X (where X is the wavelength of the radiation). 
Numerical results are presented for five different discretizations, the first of which was 
derived from a mesh (i.e., a set of patches) obtained by dividing each face into four 
squares. The second mesh was constructed from the first by subdividing each square into 
four smaller squares. The third mesh was constructed from the second by subdividing 
edge-touching patches in half along a line parallel to the edge; patches adjacent to two 
edges (i.e., corner patches) were divided into quarters. Meshes for the fourth and fifth 
discretizations were constructed by recursively applying the procedure by which the third 
mesh was constructed from the second. This process, known as patch tapering, is useful 
for resolving the source singularities that arise in the vicinity of geometric singularities. 
It also puts stress on the conventional CFIE because points near the edges get close 
together. Table 3 lists the maximum and average number of iterations the CGS routine 

required to obtain solutions for 92 independent excitations to a residual error of 1CT 3 . 
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! 

The total number of unknowns is the result of using a 9-point quadrature rule on each 
square or rectangular patch. The iteration count for CFIE2 grows very slowly with 
increasing taper depth, whereas for CFIE1 it increases steadily as expected. 


unknowns 

taper 

CPIE1 

CPIE2 



depth 

max 

ave 

max 


ave 

432 

0 

12 

6.5 

10 


4.3 

1728 

1 

18 

9.9 

11 


4.9 

3888 

2 

26 

14 

11 


4.9 

6912 

3 

41 

23 

11 


5.5 

10800 

4 

58 

36 

13 


5.7 


Table 3: Iteration count vs. taper depth for IA PEC cube. 


f: Appendix - Second-Kind Integral Operators 

B| " The standard definition of a second-kind integral operator is an operator of the form 

JSCS; 

1 10 AI + K, (38) 
fH where X is a constant, / is the identity, and iTis a compact operator. In scattering theory, 
^ one encounters operators of the form 

^+^ 2 P 2 +^, (39) 
where \ and are constants and Pi and P 2 are orthogonal projection operators such that 

(40) 

Operators of the form shown in equation (39) possess most of the desirable properties of 
second-kind integral operators. Although somewhat non-standard, such expressions are 
20 referred to as second-kind operators throughout this description. 


15 


P l+ P 2 =I. 
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Glossary of Terms and Conventions 

Composition Operator - The composition operator used herein, and represented by the 
symbol, ° , describes the following operation. Given a first integral operator, 

A s jds'K A (x 9 x f ) , and a second integral operator, B = jds"K B (x\ x") , then the 

s s 

5 composed operator AoB= fes f K A (x, x') \ds"K B (x\ x n ) . 

s s 

Equivalent surface sources - According to a well-known mathematical theorem (the 
"equivalence principle"), the fields outside a volume V can be completely specified by a 
distribution of equivalent sources on the boundary of the volume S. If S corresponds to 
the surface of a metal target (which approximates a perfect electrical conductor), the 
equivalent sources correspond to true, physical surface currents. 

Means - Generally, for purposes of this invention, the term refers to a computer program 
or sub-program, or a portion of a computer program operative for performing the task 
that the means is for. The means is generally considered to be language and computer 
platform independent, and in the form of computer instructions (or code) that may be 
either in the form of software or hardware. 

Object - This term object refers to primarily any arbitrarily shaped closed surface for 
20 which equivalent surface sources and the resultant electromagnetic scattering are to be 
determined. Objects of interest could include, as non-limiting examples, aircraft or parts 
thereof, automobiles, tanks, boats, missiles, and fixed structures. The equations set forth 
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herein are developed for perfectly electrically conducting (PEC) materials or 
approximately perfectly conducting materials such as metals. However, extension to 
surface scattering with more general boundary conditions is considered routine and 
within the scope of the present invention. The extension to an interface between two or 
more dielectrics, for example, is straightforward, with the resulting operators having 
condition numbers slightly lower than the case described herein. Although primarily 
developed for closed surfaces, the approach can potentially be extended for the case of 
open, perfectly conducting surfaces. 
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